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d . In this paper numerical methods for solving stochastic differen- 

tial equations with Markovian switching (SDEwMSs) are developed 
by pathwise approximation. The proposed family of strong predictor- 
corrector Euler-Maruyama methods is designed to overcome the prop- 
agation of errors during the simulation of an approximate path. This 
<^ paper not only shows the strong convergence of the numerical solu- 

tion to the exact solution but also reveals the order of the error under 
some conditions on the coefficient functions. A natural analogue of 
i-^H , p-stability criterion is studied. Numerical examples are given to il- 

lustrate the computational efficiency of the new predictor-corrector 
Euler-Maruyama approximation. 
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1 Introduction 



£T) \ Stochastic differential equations with Markovian switching (SDEwMSs) 

arise in mathematics models of hybrid systems that possess frequent unpre- 
dictable structural changes. One of the distinct features of such systems is 
that the underlying dynamics are subject to changes with respect to certain 
configurations. Such models have been used with great success in a vari- 
ety of application areas, including flexible manufacturing systems, electric 
, power networks, risk theory, financial engineering and insurance modeling, 

we refer the readers to Arapostathis, Ghosh and Marcus [I], Jobert and 
Rogers [7], Mao and Yuan [11], Rolski, Schmidli, Schmidt and Teugels [LI] . 
Smith [TS], Yang and Yin [IB] and references therein. 

Generally, although the fundamental theories such as existence and unique- 
ness of the solution as well as stability of SDEwMSs have been well studied, 
most of SDEwMSs cannot be solved analytically. Thus, appropriate numer- 
ical approximation methods such as the Euler (or Euler-Maruyama) method 
are needed to apply SDEwMSs in practice or to study their properties. 
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Yuan and Mao [18J firstly considered the numerical solutions of the fol- 
lowing stochastic differential equations with Markovian switching 

dy(t) = f(y(t),r(t))dt + g(y(t),r(t))dW(t), (1.1) 

here y(t) is referred to the state while r(t) is regarded as the mode. The 
system will switch from one mode to another in a random way, and the 
switching between the modes is governed by a Markov chain. They proved 
the mean-square convergence of the Euler-Maruyama(EM) approximation 
for this hybrid stochastic systems, and the order of error was also estimated. 
Yin, Song and Zhang [IT] extended (jl.ip to a family of more general jump- 
diffusions with Markovian switching, and proved the numerical solutions 
based on finite-difference procedure weak converge to the desired limit by 
means of a martingale problem formulation. 

During recent years, there also exist extensive literatures which prove 
the convergence of the Euler-Maruyama method applied to some stochastic 
differential equation with some additional feature, like including some sort of 
delay, jumps, Markovian switching or combinations thereof, see for example, 
Bruti-Liberati and Platen [2], Hou, Tong and Zhang [6j, Li and Hou [9], 
Mao and Yuan [11], Rathinasamy and Balachandran [13], among others. 
The corresponding proof is basically the same each time, the only novelty 
coming from changing it a bit to deal with the additional feature. 

It is well known that Euler-Maruyama method and most other explicit 
schemes for solving stochastic differential equations (SDEs) work unreliably 
and sometimes generate large errors, see for instance Milstein, Platen and 
Schurz [10] , implicit and predictor-corrector schemes are designed to achieve 
improved numerical stability and turn out to be better suited to simulation 
task. Generally, implicit schemes usually cost significant computational time 
and are sometimes not reliably accomplished, however, this phenomenon can 
be avoided when using some appropriate discrete time schemes, including 
predictor-corrector methods. In Kloeden and Platen [8], predictor-corrector 
methods have been proposed as weak discrete time approximations for solv- 
ing SDEs, which can be used in Monte Carlo simulation. For the strong 
discrete time approximation of solutions of SDEs, a family of predictor- 
corrector Euler methods has been developed in Bruti-Liberati and Platen 
[3]. However, there are no strong predictor-corrector methods available for 
SDEwMSs yet. 

In this paper, we develop a new family of strong predictor-corrector 
Euler-Maruyama (PCEM) methods for SDEwMSs (jl.ip . which are shown 
to converge with strong order 0.5, and demonstrate their performance by 
considering some examples. 

The rest of the paper is arranged as follows. In Section 2 we introduce 
some necessary notations and define a family of strong predictor-corrector 
Euler-Maruyama approximate solutions to SDEwMSs. In Section 3 we show 
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that the PCEM solutions converge to the exact solution in L 2 under the 
global Lipschitz condition and reveal the order of convergence is 0.5. In Sec- 
tion 4 we extend the PCEM convergence results to multi-dimensional case 
under certain conditions. In Section 5 the numerical stability of SDEwMSs 
will be introduced and discussed. Finally, in Section 6 some numerical ex- 
amples are given and compared for simulated paths with different degrees 
of implicitness to illustrate the computational efficiency of the predictor- 
corrector Euler-Maruyama approximation. 



2 Preliminary and algorithm 

Let (Q, J 7 , {J-i}t>o, -P) be a complete probability space with a filtration 
{•^i}t>o satisfying the usual conditions. Suppose that there is a finite set 
S = {1,2, ...,N}, representing the possible regimes of the environment. We 
work with a finite-time horizon [0, T] for some T > 0. Let /(■,•) : K d x S — > 
R d , g(-, ■) :M. d xS^- R dxd be both Borel measurable. Consider the dynamic 
system given by (1.1) with initial value y(0) = yo £ M. d and r(0) = iq € S, 
where W{t) = {W l {t), • • • , W d {t)) T is a d-dimensional J-^-adapted standard 
Brownian motion, and r(t) is a continuous-time Markov chain taking value 
in a finite state space S = {1,2, N} with the generator Q = (qij)NxN 
given by 



P{r(t + 5) =j\r(t) =i} = < 
provided S I 0, and 



qijS + o(S), ifi^j, 
l + q u 6 + o(S), if i = j, 



-qu = ^qij < +oo. 

i¥=3 

We assume W(t) and r(t) are independent. Throughout this paper, we 
denote by | • | the Euclidean norm for vectors or the trace norm for matrices. 

2.1 Existence and uniqueness 

Under certain conditions we can establish the existence of a pathwise 
unique solution of (jl.ip . Here we make the following global Lipschitz (GL) 
and linear growth (LG) assumptions: 

(HI) GL: For all (x, i), (y,i) € M d x S, there exists a constant L\ > such 
that 

l/OM) - f(y,i)\ 2 + -g(y,i)\ 2 <Lx\x- y\ 2 . 

(7i2) LG: For all (x, i) G W 1 x S, there exists a constant L2 > such that 
\f(x,i)\ 2 \/\g(x,i)\ 2 <L 2 {l + \x\ 2 ). 
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Remarks 2.1. It is easy to show that if /(•,•)> ff(v) satisfy GL condition, 
then they also satisfy LG condition, but for the convenience of the reader 
we preserve it. 

The following theorem guarantee the existence and uniqueness of the 
solution to equation (jl.ip . which are of use in studying some numerical 
schemes. 

Theorem 2.1. If f(x,i), g(x,i) satisfy the conditions (HI), (H2), and 
suppose W(t),r(t) be independent. Then there exists a unique d- dimensional 
J-t- adapted right- continuous process y(t) with left-hand limits which satisfies 
equation such that y(0) = yo and r(0) = iq a.s. 

Proof. See Theorem 3.13 in Mao and Yuan [11]. □ 



2.2 Algorithm 

Now we turn our attention to numerical algorithm. For convenience, 
we first consider one-dimensional SDEwMSs. Given A > as a step size, 
denote {^}j>i the usual equidistant time discretization of a bounded interval 
[0, T], i.e. to = 0, ti — = A, if t n -\ < T < t n then set t n = T. Denote 
AW tk =W(t k+1 )-W(t k ). 

For given partition {t k } k >i, {r(t k )}k>i is a discrete Markov chain with 
transition probability matrix (P(i,j))NxN, here P(i,j) = P(r(tk+i) = 
j\ r (tk) = i) is the ijth entry of the matrix e^ tk+1 ~ tk ^ , thus we could use fol- 
lowing recursion procedure to simulate the discrete Markov chain {r(tk)}k>i, 
suppose r(tk) = i\ and generate a random number £ which is uniformly dis- 
tributed in [0,1], then we define 



r(t k +i) 



i 2 , if i 2 G S - {N} and Z?=i P(h,j) < Z < E?=i P(h,3 



[N, if Eft l P(ii,j)<f. 



Repeating this procedure a trajectory of {r(tk)}k>i can be simulated. 

Now we can introduce a new family of strong predictor-corrector Euler- 
Maruyama(PCEM) methods for SDEwMSs. Given initial value Y to = y G K 
and rt = iq G 5, the proposed family of strong PCEM is given by the 
predictor 

Y tk+1 =Y tk + f(Y tk ,r tk )A + g(Y tk ,r tk )AW tk , (2.1) 

and by the corrector 

Yt k+1 =Y tk + {ef v (Y tk+1 ,r tk ) + (1 - 9)f v (Y tk ,r tk )}A 
+ iV9(Y tk+1 , r tk ) + (1 - V )g(Y tk , r tk ) } AW tk . 
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Here parameters 9, rj € [0, 1] denote the degree of implicitness in the drift 
and the diffusion coefficients, respectively, and f v (x,i) is defined as 

/„(*, i) = f(x, i) - V g(x, i) , r, € [0, 1], (2.3) 

which is called the corrected drift function. This scheme can be written in 
the form 

4 

Yt k+1 = Y tk + f(Y tk ,r tk )A + g(Y tk ,r tk )AW tk +^Ri, k , (2.4) 

1=1 

where 

Ri,k = 0{f(Y tk+1 ,r tk ) - f(Y tk ,r tk )}A, (2.5) 

R2,k = -Br]{g(Y tk+1 ,r tk ) — g{Y tk ,r tk ) — }A, W 

R 3 , k = -V9(Y tk ,r t f 9{Y ^ A, (2.7) 

Ra,u = r]{g(Y tk+1 ,r tk )- g(Y tk , r t J}A W t „ . (2.8) 
For each t € [tk,tk+i), let 

= K tfc ,f(t) = r t „,Ri(t) = ^,l = l,2,3,i? 4 (t) = (2.9) 

Therefore, we can define the continuous approximation solution Y(t) on 
the entire interval [0, T] by 



Y(t)=y + /* f(Y(s),f(s))ds+ f g{Y(s),f{s))dW{ t 
Jo Jo 
3 ft _ ft 
+ V / Ri(s)ds+ / R A {s)dW{s). 
l=l Jo Jo 



(2.10) 



Note that Y(tf.) = Y(tk) = Yf k , which means Y(t) and Y(t) coincide 
with the discrete approximate solution at the gridpoints. 

Remarks 2.2. The major advantage of the above PCEM approximate schemes 
is that there are flexible degrees of implicitness parameters 6 and ij to choose 
for simulating paths properly. For the case 9 = r] = one recovers the Euler- 
Maruyama scheme which is well discussed in Yuan and Mao [18] . 
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3 Convergence with the global Lipschitz(GL) and 
linear growth(LG) conditions 

In this section, we will prove that the numerical solution Y(t) converges 
to the exact solution y(t) in L 2 as step size A \. 0, and the order of conver- 
gence is one-half. To begin with, we need the following GL condition and 
LG condition for /?,(•, ■). 

(TiS) GL: For all (x,i), (y,i) GKxS, there exists a constant L3 > such 
that 

\frj(x,i) - fn(y,i)\ 2 < L 3 \x - y\ 2 . 
(TiA) LG: For all (x, i) € R x S, there exists a constant L4 > such that 

\f v (x,i)\ 2 <L 4 (l + \x\ 2 ). 

We are now ready to present the key results of this section which are 
stated as following. 

Theorem 3.1. Assume the SDEwMSs < TO)) defined on (O, JF , {J"t}t>o, P) 
satisfying 

a) W(t),r(t) are independent, 

b) /(•> •)>#('> -),fri(-> ■) satisfy conditions (Til), (H2), (US) and (H4), 

then the unique strong solution y(t) and numerical solution Y(t) obtained 
in section 2.2 satisfying: 

E( sup \Y(t) -y(t)\ 2 ) < CA + o(A), (3.1) 

0<t<T 

where C is a positive constant independent of A. 

In order to give the proof of this theorem, we first provide a number of 
useful lemmas. The first two lemmas show that the continuous approxima- 
tion has bounded moments in a strong sense. The latter lemmas play an 
important role in proving the strong convergence result, which mainly refer 
to Bruti-Liberati and Platen [3]. 

Lemma 3.1. Under conditions (Til), (Ti2), (TiS) and (Ti4), for any p > 
2, there exists a constant M which is dependent on p, T, L and yo, but 
independent of A, such that 

E( sup \Y(t)\ p ) < M. (3.2) 

0<t<T 

We omit the proof since it is similar to Lemma 4.1 in Mao and Yuan [llj. 
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Lemma 3.2. Under conditions {Til), ("H2) ; (*H3) and (%4), i/jere exists a 
constant M which is dependent on T, L and yo, but independent of A, such 
that 



E( max \Y t ,r) < M. (3.3) 

0<fc<n T k 

This is an immediate result of Lemma 13. 11 since Y(tf~) = Y(t)~) = Yt k . 

Lemma 3.3. There exists a constant C which is dependent on T, L and yo, 
but independent of A, such that 

4 

V£[ max | V R hk \ 2 ] < CA. (3.4) 

* — ' Kn<riT — ' 
1=1 ~ ~ 0<fc<n-l 

Proof. By the Cauchy-Schwarz inequality and the GL condition, from equa- 
tion (|2,5p . we can obtain 

|2l 

i<riT 

0<fc<n-l 



£/[ max j N 

Kn<nT ^ — ^ 
~ ~ 0<fc<n-l 

= E[ m f | £ 0{/(y, fc+1 ,r,J-/(y tfc ,r t J}A| 2 ] 

~ ~ 0<fc<n-l 

<M max ( l 0A l 2 )( E -m h ,rt k )\ 2 )] 

Kn<jiT ' — ' ' — ' 



i<riT 

0<fe<n-l 0<fc<n-l 

<CAE[{ A )( E \f& k+1 ,rt k )-f(Y tk ,r tk )\ 2 )} 

0<fc<n T -l 0<fc<n T -l 

<CA£[ £ \Y tk+1 -Y tk \ 2 ] 

0<fc<n T -l 

<CA£[ ^ (£;[|/(F tfe ,r tfc )A| 2 |^J+^[| 3 (F t „ nfe )A^J 2 |J- t J)]. 

0<fc<n T _! 

(3.5) 

Then by using the Cauchy-Schwarz inequality, the Ito's isometry, the LG 
condition and Lemma 13.21 we get 

E[ max j Rik\ 2 ] 

KtI<IIT ^— 

~ ~ 0<fc<n-l 



<CAE[[ tnT E((l + \Y tn f)\T tnz )dz] 
Jo 

rtn T 

< CA / E(l+ max \Y tn \ 2 )dz 

Jo 0<n<n T 



(3.6) 



< CA. 

With the similar steps as in (|3.5p and (|3.6p . we have 
£7[ max | V i? 2 fc! 2 ] < CA. 

Kn<riT * — ' 



1<1%t 

0<k<n-l 



(3.7) 
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It is also easy to show by the LG condition and Lemma 13, 21 that 

^fel E ^, fc | 2 ]< CA. (38) 

~ ~ 0<fc<n-l 

By using Doob's martingale inequality, the Ito's isometry, the GL condition 
and equation (|2.8p . we have 



t fe +i 



E\ max i > Rd k\ 

~ ~ 0<fc<rt-l 

= E[ max | ^ W" ( 5 (F tfc+1 ,r t J- ff (y t)b ,r t J)^)p 

Kn<riT ' — ' /j- 
~ ~ 0<fc<n-l Jtk 

< CE[\ [ T (g(Y tnz+1 ,r tn ) - g(Y tnz ,r tn ))dW(z)\ 2 } ( 3 - 9 ) 



o 



T 



C / S[| 5 (F w+1) r w )- 5 (F w ,r w )|1^ 



(3.10) 



<c/ s[|y tnz+1 -y tn j 2 ]dz. 

Jo 

Therefore, with similar steps as in (|3.5p and (|3.6p . we also have 
E\ max I > i?A i, | 2 1 

L Kn<n T ' ^ ' 1 1 

- - 0<fc<n-l 

<C f E[ f nZ+1 E[(l + \Y tnz \ 2 )\T t Jds]dz 
Jo Jt nz 

<C [ T E[[ tnZ+1 E[(l+ max \Y tn \ 2 )\T tn Jds]dz 

Jo Jt nz 0<n<n T 
< CA. 

Thus the required assertion follows. □ 

Lemma 3.4. Under conditions (HI), (H2), (713) and (HA), then for any 
t £ [tk-,tk+i), we have 

J>[ sup I [ S ^dz\ 2 ]+E[ sup | [ S ^dW(z)\ 2 ]<o(A). (3.11) 

^ t k <s<t Jt k ^ t k <s<t Jt k ^W tk 

The proof of this lemma is similar to that in Lemma 13.31 
Now we are in a position to prove our Theorem 13. 11 
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Proof of Theorem 13.11 : From equation (|2,10p , we have 
Z(t) = E( sup \Y( S )-y(s)\ 2 ) 

0<s<t 

= E[su V | [ S (f(Y(z),r(z))-f(y(z),r(z)))dz 
0<s<t Jo 

+ f S (g(Y(z),f(z)) - g(y(z),r(z)))dW(z) ( ' > " l2) 
Jo 

+ V / Ri(z)dz+ / R 4 (z)dW(z)\ 2 }. 
l=1 Jo Jo 

Let n = [t/A], the integer part of t/A. Then, by Holder inequality, Doob's 
martingale inequality and equation (|2.9p . we have 

Z(t)=E( sup \Y(s)-y(s)\ 2 ) 

0<s<t 

<CE f \f(Y(z),f(z)) - f(y(z),r(z))\ 2 dz 
t 

2. 







+ CE / \g(Y(z),r(z))-g(y(z),r(z))\ 2 dz 
Jo 

4 3 



+ max | V R lk \ 2 ] + C sup | / 

^ — ' Km<n — ' * — ' < Q <+ /j- 

+ sup | f^dW(z)\ 2 ]. 

t k <s<t Jt k ^ W t k 



^dz\ 2 ] 



(3.13) 



We focus on the last three parts of the right side. From Lemma 13.31 and 
Lemma 13^41 we have 



4 3 . 

C^Mmax | £ it^| 2 ] + C^£[ sup | / 

1 — * Km<n * — * * — ' +, /+, 

1=1 ~ ~ 0<k<m-l 1=1 l k^s^i Jt, 

+ CE[ sup | [ S ^fdW{z)\ 2 ]<CA + o(A). 



SRl > k dz\ 2 ] 



A 



(3.14) 



Let Iq be the indicator function for set G. Let t € [tk,tk+ij- From (|2. 1 j) 
(|2.10p . obviously, we have Y(ty.) and I{r(Mr(t k )\ are conditionally indepen- 
dent with respect to the cr-algebra generated by r(tfc). So by the same 
procedures as in Theorem 3.1 in Yuan and Mao |18| . we can obtain 



t 

2, 



E / \f(Y(s),f(s))-f(y(s),r( S ))\ 2 ds 

J ° t (3-15) 

<2L 2 E\Y(s)-y(s)\ 2 ds + CA + o(A), 
Jo 



9 



E [ \b(Y(s),f( S )) -b(y(s),r(s))\ 2 ds 
Jo 



(3.16) 



<2L 2 [ E\Y(s)-y(s)\ 2 ds + CA + o(A). 
Jo 

Substituting ([513]) . (ETT51) . (I3TT6]) into (|3"7T3l shows that 

£?( sup |y(s) -y(s)| 2 ) < C / £?|?(a) - y(s)\ 2 ds + CA + o(A). (3.17) 

0<s<t JO 

Note that 

E\Y(s) - y(s)\ 2 < 2E\Y{s) - y(s)\ 2 + 2E\Y(s) - Y(s)\ 2 . (3.18) 

Suppose tk < s < ffc+i, by (|2.10p . Lemma [331 and Lemma |3~H we can then 
show in the same way as in the case of SDEs that 

E\Y(s) - Y(s)\ 2 < CA. (3.19) 

Substituting (|3.18p . (|3.19j) . into ()3.17|) immediately shows that 

E{ sup \Y{s)-y{s)\ 2 ) 

0<s<t 

< C J\e\Y{s) ~ V(s)\ 2 + E\Y(s) - Y(s)\ 2 )ds + CA + o(A) (3 . 20 ) 
<C I E(sup \Y(s)-y(s)\ 2 )ds + CA + o(A). 

Jo 0<s<t 

Therefore, from Gronwall inequality we obtain that 

E{ sup \Y(t)-y{t)\ 2 ) <CA + o(A). 

0<t<T 

The proof is complete. 

4 The general multi-dimensional case 

The results derived in Section 3 can be easily generalized to the multi- 
dimensional case, we just summarize the related numerical schemes and the 
convergence results in this section. 

Consider the solution y(t) = {(y 1 (t), y d (t)) T , t > 0} of the d-dimensional 
SDEwMSs 

ft m ft 

y(t) = y(P)+ / f(y(s),r(s))ds + J2 / gi(y(s),r(s))dW3(s), 
Jo j=1 Jo 

for t > 0. Here y(0) € M. d denotes the deterministic initial value, = 
{W 3 (t),t > 0}, j G {1,2, ...,m} is a standard Brownian motion. r(t) is a 



10 



Markov chain. The function / : M. d x S h4 M d has the fcth component •). 
The function ^ : M. d xS \— > M. d , j 6 1, 2, m has the fcth component g '■*(■, •). 
We define the function /„ for 77 € [0, 1] with the fcth component 

771 d r) "i 32 ( '\ 

fH{x,i) = f k {x t i)-7 fk Y.»'- J '■->■■< ' ' d J - 

jiJ2=l i=l 

for (x, i) € x 5, which satisfies the following GL condition and LG con- 
dition 

(H3') GL: For all (x, i), (y, i) € M. d x 5, there exists a constant L3 > such 
that 

|/»?0M) - fr,(y,i)\ 2 < L 3 \x - y\ 2 . 
(H4f) LG: For all (x, i) € M. d x S, there exists a constant L4 > such that 

|/„(M)I 2 <l 4 (i + |x| 2 ). 

The kth. component of the proposed family of strong PCEM schemes is 
given by the predictor 

m 

Y£ +1 = Yt + f h (Yt n ,r tn )A + Y,9 k ' j (Y tn ,r tn )AWi, 

and by the corrector 

Y t k n+1 =Yt + {e k fi(Y tn+1 ,r tn ) + (1 - e k )f^Y tn , rtn )}A 

m 

+ J2^ j (Y tn+1 , rtn ) + (1 - ^)^'(y tn) r tn )}AW tn) 
3=1 

for ^,r/ fc G [0,1], € {1, 2, d}. 

Hence we can define the approximation solution Y(t) similarly to equa- 
tion (|2,1U|) . then we can derive the following theorem analogically. 

Theorem 4.1. Assume the SDEwMSs defined on (Q, T, {Tt}t>o, P) 

satisfying 

a) W(t),r(t) are independent, 

b) /(•,•)> 9(-,-),f v (-r) satisfy conditions (HI), {H2), (H5t ') and (Hi), 

then the unique strong solution y(t) and numerical solution Y(t) satisfying: 
E{ sup \Y(t)-y{t)\ 2 ) <CA + o(A), 

0<t<T 

where C is a positive constant independent of A. 
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5 Numerical Stability 



In this section we consider numerical stability issues, extending the anal- 
ysis in Platen and Shi [12J to the Markovian switching case. When simulating 
discrete time approximations of solutions of SDEwMSs, numerical stability 
is clearly as important as numerical efficiency There have been various ef- 
forts made in the literature trying to study numerical stability for a given 
scheme approximating solutions of SDEs, see, for instance, Hofmann and 
Platen [5], Higham [JJ, Bruti-Liberati and Platen [3] and Platen and Shi 
|12| . Generally for analyzing numerical stability, some specifically designed 
test equations are necessary to be introduced, the test SDEs used in the 
above literatures are linear SDEs with multiplicative noise defined as 

dX t = (1 - ^a)XX t dt + y^\X\X t dW t , (5.1) 

for every t > 0, where Xq > 0, A < and a G [0, 1). Its explicit solution is 
of the form 

X t = X exp{(l - a)Xt + ^/a\X\W t }, t > (5.2) 

As a unified criterion, Platen and Shi [12] proposed the concept of p- 
stability criterion, which means that a process is p-stable if in the long run 
its pth moment vanishes. Hence, for p > 0, A < 0, p-stable in equation (|5.ip 
may be characterized by 

lim E{\X t \ p ) = iff < a < — 

t-i-OO 1 + ^ 

Since for different combinations of values of A A, a and p with given time 
step size A, a discrete time approximation and the original continuous 
process X t have different stability properties. To explore these differences, 
the concept of stability region is introduced. The stability region, denoted 
by T, is determined by those triplets (XA,a,p) G (— oo,0) x [0, 1) x (0, oo) 
for which the discrete time approximation Yt is p-stable with time step size 
A, when applied to the test equation ()5.ip . 
By defining the random variable 

G n+1 (XA,a) = [-^|, 

* n 

which is called the transfer function of the approximation Yt at time t n , 
Platen and Shi [12j derive the following useful result which can determine 
the stability regions for given schemes by the following theorem. 

Theorem 5.1. A discrete time approximation is for given X < 0, a G [0, 1) 
and p > 0, p-stable if and only if 

E((G n+l (XA,a)r) < 1. 
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In the spirit of Platen and Shi [12J, stability regions for a range of schemes 
of SDEwMSs are discussed in this paper. Here we consider the test process 
X r (t) = {Xt, r (t) ) t > 0} satisfies the linear SDEwMSs with multiplicative 
noise 

dX t = (1 - ~a(r(t)))X(r(t))X t dt + y/a(r(t))\X(r(t))\X t dW t , (5.3) 

for every t > 0, where r(t) is a Markov chain taking values in a finite state 
space S = {1, 2, JV}, r(0) =ioeS,Io = * r (o) > and (a(r(t)), A(r(t))) G 
{(a(i),A(i)),a(i) G [0,1), A(i) < 0,i = 1,2, ...,JV}. It is well known that the 
explicit solution of the test equation (|5.3p is (see Mao and Yuan [11] ) 

X t = X exp{ [ (1 - a(r(t)))A(r(t))£ft + f y/a(r(t))\X(r(t))\dW t }. (5.4) 
J o J o 

For the convenience of numerical comparison, we introduce the following 
stability criterion: 

Definition 5.1. For p > 0, a process Y r n\ = {Y t>r M,t > 0} is called state- 
p-stable if for each i = 1, 2, N, Yj, = {Y tt i, t > 0} satisfies 

hm ^(|y M | p ) = 0. 

t— >oo 

Definition 5.2. The state- stability region T s is determined by those triplets 
(AA, a,p) € (— oo, 0) x [0, 1) x (0, oo) for which the discrete time approxima- 
tion Y r M is state-p-stable, when applied to the test equation A5.3\) . for each 
i = 1,2, ...,N, (ct(i), X(i)A,p) G T s with time step size A. 

Then we can obtain the following conclusion: 

Theorem 5.2. The state-stability region T s , which is generated by one algo- 
rithm applied to the test equation 15. 3\) . is the same as the stability region T, 
which is generated by the same algorithm applied to the test equation \5.1\) . 

Proof, r C T s is obvious. To prove T s C T, suppose for all i = 1,2, ...,iV, 
(a(i), \(i)A,p) G T s . By Definition 15.21 we can see that Y r n) is state-p- 
stable. And by Definition 15 . 1 1 we know for each i = 1, 2, N,Yi = {Ytj, t > 
0} = Pt,(a(i),A(i)A,p)j* > 0} is P-Stable, so (a(i), X(i)A,p) G T, for all i = 
1,2, ...,N. Immediately, we have r s C T. □ 

Remarks 5.1. By the conclusions derived in Platen and Shi [12], we can also 
see that the PCEM methods are more efficient than the EM method under 
these conditions in SDEwMSs case. 
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6 Numerical examples 



In this section, we discuss two numerical examples to illustrate our theory 
established in the previous sections. Let us now consider several combina- 
tions of parameters and rj in equation ()2.2|) . the names of the methods 
listed below are similar to those used in Bruti-Liberati and Platen [3] and 
Platen and Shi [12]. 

(1) = 0, 7] = (called EM scheme), 

(2) = hiV = \ (called symmetric PCEM scheme), 

(3) =2)^ = (called semi-drift-implicit PCEM scheme), 

(4) = 1, 77 = (called drift-implicit PCEM method), 

(5) = 0, rj = § (called semi-diffusion-implicit PCEM scheme), 

(6) = 1, 77 = 1 (called fully implicit PCEM scheme). 

For a given problem, we will compare the simulated paths for these 
different degrees of implicitness. If these paths differ significantly from each 
other, then some numerical stability problem is likely to be present and one 
needs to make an effort in providing extra numerical stability for further 
research. 

Example 6.1. Let W(t) be a scalar Brownian motion. Let r(t) be a right- 
continuous Markov chain taking values in S = {1,2} with generator 



Q 



-1 1 

q -q 



And W(t) and r(t) are assumed to be independent. Consider an one- dimensional 
linear SDWwMS 

dy(t) = y(t)a(r(t))ds + y(t)b(r(t))dW(t) (6.1) 

for t > 0, where 

o(l) = 1, 6(1) = 2 
o(2) = 2, 6(2) = 1 



It is well known that equation (|6.1|) has an explicit solution 



Jo 2 j 



y{t) = y exp[ / a{r(s))ds + I (r{s))dW{s) - £ / b\r(s))ds}. (6.2) 



For simulation reason, it is convenient to transform (|6.2p into following 
recursion form with y(ta) = yo, 



y{tk+i) ~ y(t k )exp[(t k+1 - t k )a{r tk ) + (W(t k+1 ) - W(t k ))b(r tk 
1 

~ 2 



1 ( 6 - 3 ) 

-{t k+1 -t k )b\r tk )\- 
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Notice that y{t k+ i) in ()6,3p is not the exact value of y(t) at the division 
points ifc+i, because r(s) is not necessarily constant on [t k ,t k +i]- However, 
since 

P{r(i fe+1 ) = i|r(tfe) =i} = l + qu(t k+ i - t k ) + o(t k +i - t h ) ->• 1 

as A 4 0, for sufficiently small A, it is reasonable to use (|6.3p as an approx- 
imation of the exact solution of y(t). 

Case 1. Let q = 2, y = 200, r = 1, A = 10~ 5 , fc = 0,1,..., 5 x 10 6 , 
namely for the corresponding time < t < 50. Compute the one-step 
transition probability matrix 

- [°- 99999 0.00001" 
^ ' ~ [0.00002 0.99998 

for the discrete Markov chain rt h = r(fcA). 

By applying the previously described procedure, the trajectory of the 
approximate solution Y(t) with given stepsize A can be constructed. In 
this paper, we do not draw the figure of the simulating trajectory, instead, 
to carry out the numerical simulation clearly, we repeatedly simulate and 

compute sup (\Yrg. Vi )(tk) —y{tk)\ 2 ) (i = 1, 2, 3, 4, 5, 6) for 1000 times, then 

t fc e[o,50] 

calculate the sample mean E( sup \Y(g. Vi }(tk) —y{tk)\ 2 )(i = 1,2,3,4,5,6). 

i fc e[0,50] 

The results are listed in the following Table 1. 



di,Vi 


E( sup \Y { g it1H) (t k ) -y(t k )\ 2 ) 

t fc G[0,50] 


0,0 


1.90674403017316e+30 


1 1 

2' 2 


7.84001334096008e+25 


i,0 


1.90607879136196e+30 


1,0 


1.90541425316912e+30 


0,i 


7.21187503884961e+25 


1,1 


1.99576630102430e+30 



Table 1. Estimation of the errors between the numerical solutions and exact solution 



Case 2. Let q = 1.5, y = 200, r = 1, A = 10" 5 , k = 0,1,..., 5 x 10 6 , 
namely for the corresponding time < t < 50. Compute the one-step 
transition probability matrix 

(A _ I" 0.99999 0.00001 " 
^ ' ~ [0.000015 0.999985 

for the discrete Markov chain rt k = r(kA). 

To carry out the numerical simulation we repeatedly simulate and com- 
pute sup (\Y( d . iV .)(tk) - y{t k )\ 2 ) (i = 1,2,3,4,5,6) for 1000 times, then 

t fc G[0,50] 
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calculate the sample mean E( sup \Y^g. Vi ^(t k ) —y(t k )\ 2 )(i = 1,2,3,4,5,6). 

* fc e[0,50] 

The results are listed in the following Table 2. 





E( sup \Y^ g .^(t k ) - y(t k )\ 2 ) 
t k e[o,50] 


0,0 


3.33550923616514e+36 


i i 

2' 2 


1.71565111220697e+33 


i,0 


3.34748775539284e+36 


1,0 


3.35948703848527e+36 


0.3 


3.51753180196352e+31 


1,1 


3.23681904239271e+36 



Table 2. Estimation of the errors between the numerical solutions and exact solution 



Case 3. Let q = 0.5, y = 200, r = 1, A = 10" 5 , k = 0, 1, 5 x 10 6 , 
namely for the corresponding time < t < 50. Compute the one-step 
transition probability matrix 

_ I" 0.99999 0.00001 " 
^ ' ~ [0.000005 0.999995 

for the discrete Markov chain rt k = r(kA). 

To carry out the numerical simulation we repeatedly simulate and com- 
pute sup {\Y/g um )(t k ) - y(t k )\ 2 ) (i = 1,2,3,4,5,6) for 1000 times, then 

t fc e[o,50] 

calculate the sample mean E( sup \Yrg. Vi )(tk) — y{t k )\ 2 ){i = 1,2,3,4,5,6). 

* fc e[o,50] 

The results are listed in the following Table 3. 



0i,Vi 


E( sup |y (e . )J? .)(t fe ) - y{t k )\ 2 ) 

t fc G[0,50] 


0,0 


1.44025299136110e+71 


1 1 

2' 2 


4.66862461497885e+66 


2>0 


1.39047511677566e+71 


1,0 


1.34154643977816e+71 


0.5 


5.41110836832358e+67 


1,1 


1.42636805044311e+71 



Table 3. Estimation of the errors between the numerical solutions and exact solution 



This example has been studied in Yuan and Mao [18] in which Case 1, 
Case 2 and Case 3 represent three different exponential stability or insta- 
bility situations respectively. However, since we can easily show that the 
equation (|6.ip does not satisfy the conditions of equation (|5.3|) . it will be 
not state-p-stable, the computer simulation results in Case 1, Case 2 and 
Case 3 illustrate this point in some extent. Nevertheless, when the parame- 
ters 9 and r] are selected reasonably, some PCEM methods are much more 
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efficient than the EM method in a certain extent. 



Example 6.2. Let r(t) be a right- continuous Markov chain taking values in 
S = {1, 2} with generator 



Q 



-0.5 
0.5 



0.5 
-0.5 



QA 



therefore, the one step transition probability from r(t k ) to r{t k+ i) is e 
Consider the same 1-dimensional linear SDWwMS 

dy(t) = y{t)a{r{t))ds + y(t)b(r{t))dW{t) (6.4) 
on t > 0, this time let a(r(t)), 6(r(i))take values as follows 

a(l) = 0.15, 6(1) = 0.1 

a(2) = 0.05, 6(2) = 0.1 

Choose initial values yo = 10, ro = 1, T is fixed at 10. By applying the 
previously described procedure, the trajectory of the approximate solution 
Y(t) with given stepsize A can be constructed. 

To carry out the numerical simulation we successively choose the stepsize 
A as the following Table 4, and for each A, we repeatedly simulate and com- 
pute sup (|Y(0. i7? .)(4) -y(t k )\ 2 ), (i = 1,2,3,4,5,6) for 1000 times, then 
t h e[o,io] 

calculate the sample mean E{ sup \Y^. Vi )(tk) — 2/(^fc)| 2 )(^ = 1) 2, 3, 4, 5, 6). 

t fc 6[0,10] 

The results are listed in the following Table 4. 



E( sup \Y i8i)Vt) (t k ) - y(t k )\ 2 ) 

t fc 6[0,10] 


A \ 9i,in 


0,0 


i i 

2' 2 


i° 


1,0 


0,i 


1,1 


0.1 


9702.4683 


15.5672 


1421.4360 


8805.6057 


7629.3505 


7912.6245 


0.01 


376.9640 


0.2299 


265.3262 


457.8189 


150.4551 


342.0570 


0.001 


24.0510 


0.0017 


23.6959 


25.8235 


1.1710 


23.8747 


0.0001 


3.0465 


2.16e-05 


3.1024 


3.1944 


0.0158 


3.0511 


0.00001 


0.1968 


1.48e-07 


0.1969 


0.1971 


9.73e-05 


0.1968 



Table 4. Estimation of the errors between the numerical and exact solutions 



Clearly, we can easily show that the equation (|6.4|) satisfies the conditions 
of equation (|5.3p . so it will be state-p-stable, and the simulation results listed 
in Table 4 just illustrate the stability properties in a certain extent. On the 
one hand, the numerical method reveals that the numerical solution Y(t) 
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defined by the strong PCEM methods converge to the exact solution y(t) in 
L 2 as step size A | 0, and the order of convergence is one-half, i.e. 



E{ sup \Y(t)-y(t)\ 2 ) <CA + o(A). 

0<t<T 

On the other hand, when the parameters selection are reasonable, the PCEM 
methods is much more efficient than the EM method, which strongly demon- 
strate our results. 
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